function graphs_data_sim(fp,col_long,data,sim_data, paper_dir)

cd(paper_dir)

% Calculate some averages by child age
xvec = [3:16]'; % kid ages
yvec = NaN*ones(length(xvec),7);
for i = 1:size(xvec,1)
    age_i = xvec(i);
    
    % 1. Parental consumption c, sim
    rows = sim_data(:,col_long.age)==age_i;
    yvec(i,1) = nanmean(sim_data(rows,col_long.c));
    
    % 2. Child investment e, sim
    rows = sim_data(:,col_long.age)==age_i;
    yvec(i,2) = nanmean(sim_data(rows,col_long.e));
    
    % 3. Child consumption x, sim
    rows = sim_data(:,col_long.age)==age_i;
    yvec(i,3) = nanmean(sim_data(rows,col_long.x));
    
    % 4. Total household income, sim
    rows = sim_data(:,col_long.age)==age_i;
    yvec(i,4) = nanmean(sim_data(rows,col_long.Y));
    
    % 5. Leisure of household members, sim
    rows = sim_data(:,col_long.age)==age_i;
    yvec(i,5) = nanmean(sim_data(rows,col_long.l1)); % mom
    yvec(i,6) = nanmean(sim_data(rows,col_long.l2)); % dad
    yvec(i,7) = nanmean(sim_data(rows,col_long.lc)); % child
end

% Alternatively: aggregate it by slightly wider age categories of 2 years to make smoother plots
xvec2_l = [3:2:15]'; % kid ages, lower bound
xvec2_h = [4:2:16]'; % kid ages, upper bound
yvec2 = NaN*ones(length(xvec2_l),16);

for i = 1:size(xvec2_l,1)
    age_lb = xvec2_l(i);
    age_ub = xvec2_h(i);
    
    % 1.a) mom hours (when working), data
    rows = (data(:,col_long.age)>=age_lb) & (data(:,col_long.age)<=age_ub) & (data(:,col_long.mom_hours)>0);
    yvec2(i,1) = nanmean(data(rows,col_long.mom_hours));
    % 1.b) mom hours (when working), sim
    rows = (sim_data(:,col_long.age)>=age_lb) & (sim_data(:,col_long.age)<=age_ub) & (sim_data(:,col_long.mom_hours)>0);
    yvec2(i,2) = nanmean(sim_data(rows,col_long.mom_hours));
    % 2.a) dad hours (when working), data
    rows = (data(:,col_long.age)>=age_lb) & (data(:,col_long.age)<=age_ub) & (data(:,col_long.dad_hours)>0);
    yvec2(i,3) = nanmean(data(rows,col_long.dad_hours));
    % 2.b) dad hours (when working), sim
    rows = (sim_data(:,col_long.age)>=age_lb) & (sim_data(:,col_long.age)<=age_ub) & (sim_data(:,col_long.dad_hours)>0);
    yvec2(i,4) = nanmean(sim_data(rows,col_long.dad_hours));
    
    % 3.a) mom proportion employed, data
    temp1 = nansum(data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub & data(:,col_long.mom_hours) > 0);
    temp2 = nansum(data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub & data(:,col_long.mom_hours) == 0);
    yvec2(i,5) = temp1/(temp1 + temp2);
    % 3.b) mom proportion employed, sim
    temp1 = nansum(sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub & sim_data(:,col_long.mom_hours) > 0);
    temp2 = nansum(sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub & sim_data(:,col_long.mom_hours) == 0);
    yvec2(i,6) = temp1/(temp1 + temp2);
    % 4.a) dad proportion employed, data
    temp1 = nansum(data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub & data(:,col_long.dad_hours) > 0);
    temp2 = nansum(data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub & data(:,col_long.dad_hours) == 0);
    yvec2(i,7) = temp1/(temp1 + temp2);
    % 4.b) dad proportion employed, sim
    temp1 = nansum(sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub & sim_data(:,col_long.dad_hours) > 0);
    temp2 = nansum(sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub & sim_data(:,col_long.dad_hours) == 0);
    yvec2(i,8) = temp1/(temp1 + temp2);
%     % 5.a) proportion both employed, data
%     temp1 = nansum(data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub & data(:,col_long.mom_hours) > 0 & data(:,col_long.dad_hours) > 0);
%     temp2 = nansum(data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub & (data(:,col_long.mom_hours) == 0 | data(:,col_long.dad_hours) == 0));
%     yvec2(i,9) = temp1/(temp1 + temp2);
%     % 5.b) proportion both employed, sim
%     temp1 = nansum(sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub & sim_data(:,col_long.mom_hours) > 0 & sim_data(:,col_long.dad_hours) > 0);
%     temp2 = nansum(sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub & (sim_data(:,col_long.mom_hours) == 0 | sim_data(:,col_long.dad_hours) == 0));
%     yvec2(i,10) = temp1/(temp1 + temp2);
    %6.a) mom active time with child, data
    rows = data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub;
    yvec2(i,11) = nanmean(data(rows,col_long.tau1));
    %6.b) mom active time with child, sim
    rows = sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub;
    yvec2(i,12) = nanmean(sim_data(rows,col_long.tau1));
    %7.a) dad active time with child, data
    rows = data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub;
    yvec2(i,13) = nanmean(data(rows,col_long.tau2));
    %7.b) dad active time with child, sim
    rows = sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub;
    yvec2(i,14) = nanmean(sim_data(rows,col_long.tau2));
    %8.a) child productive alone time with child, data
    rows = data(:,col_long.age)>=age_lb & data(:,col_long.age)<=age_ub;
    yvec2(i,15) = nanmean(data(rows,col_long.tau_c));
    %8.b) child productive alone time with child, sim
    rows = sim_data(:,col_long.age)>=age_lb & sim_data(:,col_long.age)<=age_ub;
    yvec2(i,16) = nanmean(sim_data(rows,col_long.tau_c));
end

%% Time investments by child age

% Mom active time
f = figure;
plot(xvec2_l, yvec2(:,11),'k-',xvec2_l, yvec2(:,12),'k--');
set(f,'visible','off');
xlabel('Child''s Age');
ylabel('Hours per week');
axis([2 16 0 40]);
set(gca,'XTick',xvec2_l)
h = legend('Mother''s Active Time $\tau_1$ - Data','Mother''s Active Time $\tau_1$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_tau1_smooth', 'png')
saveas(f,'Fig_tau1_smooth', 'eps')

% Dad active time
f = figure;
plot(xvec2_l, yvec2(:,13),'k-',xvec2_l, yvec2(:,14),'k--');
set(f,'visible','off');
xlabel('Child''s Age');
ylabel('Hours per week');
axis([2 16 0 40]);
set(gca,'XTick',xvec2_l)
h = legend('Father''s Active Time $\tau_2$ - Data','Father''s Active Time $\tau_2$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_tau2_smooth', 'png')
saveas(f,'Fig_tau2_smooth', 'eps')

% Child productive alone time
f = figure;
plot(xvec2_l, yvec2(:,15),'k-',xvec2_l, yvec2(:,16),'k--');
set(f,'visible','off');
xlabel('Child''s Age');
ylabel('Hours per week');
axis([2 16 0 15]);
set(gca,'XTick',xvec2_l)
h = legend('Child''s Self-Investment Time $\tau_c$ - Data','Child''s Self-Investment Time $\tau_c$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_tauc_smooth', 'png')
saveas(f,'Fig_tauc_smooth', 'eps')

%% Parental labor Supply - Intensive and extensive margin

% Mom hours when working
f = figure;
plot(xvec2_l, yvec2(:,1),'k-',xvec2_l, yvec2(:,2),'k--');
set(f,'visible','off');
% title('Mom labor supply as a function of child age (if working), data vs. sim')
xlabel('Child''s Age');
ylabel('Hours per week');
axis([2 16 0 60]);
set(gca,'XTick',xvec2_l)
h = legend('Mother''s Labor Hours $h_1$ - Data','Mother''s Labor Hours $h_1$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_h1_smooth', 'png')
saveas(f,'Fig_h1_smooth', 'eps')

% Dad hours when working
f = figure;
plot(xvec2_l, yvec2(:,3),'k-',xvec2_l, yvec2(:,4),'k--');
set(f,'visible','off');
xlabel('Child''s Age');
ylabel('Hours per week');
axis([2 16 0 60]);
set(gca,'XTick',xvec2_l)
h = legend('Father''s Labor Hours $h_2$ - Data','Father''s Labor Hours $h_2$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_h2_smooth', 'png')
saveas(f,'Fig_h2_smooth', 'eps')

% Mom proportion employed
f = figure;
plot(xvec2_l, yvec2(:,5),'k-',xvec2_l, yvec2(:,6),'k--');
set(f,'visible','off');
xlabel('Child''s Age');
ylabel('Proportion of Mothers Working');
axis([2 16 0 1]);
set(gca,'XTick',xvec2_l)
h = legend('Mother''s LFP $Pr(h_1>0)$ - Data','Mother''s LFP $Pr(h_1>0)$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_LFP1_smooth', 'png')
saveas(f,'Fig_LFP1_smooth', 'eps')

% Dad proportion employed
f = figure;
plot(xvec2_l, yvec2(:,7),'k-',xvec2_l, yvec2(:,8),'k--');
set(f,'visible','off');
xlabel('Child''s Age');
ylabel('Proportion of Fathers Working');
axis([2 16 0 1]);
set(gca,'XTick',xvec2_l)
h = legend('Father''s LFP $Pr(h_2>0)$ - Data','Father''s LFP $Pr(h_2>0)$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_LFP2_smooth', 'png')
saveas(f,'Fig_LFP2_smooth', 'eps')

%% Simulated Expenditures and Leisure by child age

% Consumption and expenditures
f = figure;
plot(xvec, yvec(:,1),'k-', xvec, yvec(:,2),'k--', xvec, yvec(:,3),'k-.',xvec, yvec(:,4),'k-o')
set(f,'visible','off');
xlabel('Child''s Age');
ylabel('Dollars per week');
h = legend('Parental consumption $c$ - Simulated','Child expenditures $e$  - Simulated', 'Child consumption $x$  - Simulated', 'Total income $Y$  - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_cexY', 'png')
saveas(f,'Fig_cexY', 'eps')

% Leisure of all members, simulated only
f = figure;
plot(xvec, yvec(:,5),'k-',xvec, yvec(:,6),'k--',xvec, yvec(:,7),'k-.')
set(f,'visible','off');
xlabel('Child age');
ylabel('Hours per week');
xlim([0 17]);
h = legend('Mother leisure $l_1$ (sim)','Father leisure $l_2$ (sim)','Child leisure $l_c$ (sim)','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_leisures', 'png')
saveas(f,'Fig_leisures', 'eps')

%% Simulated CCT use by child age parental education group
% Plot CCT, reward elasticity r, intercept b, consumption x

xvec = [3:16]'; % kid ages
yvec_bins = zeros*ones(length(xvec),3); % save CCT results
zvec_bins = zeros*ones(length(xvec),3); % save r results (conditional on CCT = 1)
wvec_bins = zeros*ones(length(xvec),3); % save b results (unconditional)
vvec_bins = zeros*ones(length(xvec),3); % save x results (unconditional)
tmp_educbins = [12 ; 16];

for i = 1:size(xvec,1)
    age_i = xvec(i);
    
    % select rows for age i, educ bin 1/2/3
    rows_educbin1 = (sim_data(:,col_long.age)==age_i) & (sim_data(:,col_long.dad_educ) <= tmp_educbins(1)); % 12 or less
    rows_educbin2 = (sim_data(:,col_long.age)==age_i) & (sim_data(:,col_long.dad_educ) > tmp_educbins(1)) & (sim_data(:,col_long.dad_educ) <= tmp_educbins(2)); % 13,14, 15 or 16 (some college or college)
    rows_educbin3 = (sim_data(:,col_long.age)==age_i) & (sim_data(:,col_long.dad_educ) > tmp_educbins(2)); % 17
    
    % educ bin 1
    temp1 = nansum(sim_data(rows_educbin1,col_long.CCT) == 1);
    temp0 = nansum(sim_data(rows_educbin1,col_long.CCT) == 0);
    yvec_bins(i,1) = temp1/(temp1 + temp0);
    
    temp1 = (rows_educbin1 == 1) & (sim_data(:,col_long.CCT) == 1);
    zvec_bins(i,1) = nanmean(sim_data(temp1,col_long.r));
    
    temp1 = (rows_educbin1 == 1);
    wvec_bins(i,1) = nanmean(sim_data(temp1,col_long.b)); % unconditional on CCT use
    vvec_bins(i,1) = nanmean(sim_data(temp1,col_long.x)); % unconditional on CCT use
    
    % educ bin 2
    temp1 = nansum(sim_data(rows_educbin2,col_long.CCT) == 1);
    temp0 = nansum(sim_data(rows_educbin2,col_long.CCT) == 0);
    yvec_bins(i,2) = temp1/(temp1 + temp0);
    
    temp1 = (rows_educbin2 == 1) & (sim_data(:,col_long.CCT) == 1);
    temp2 = sim_data(temp1,col_long.r);
    zvec_bins(i,2) = nanmean(temp2);
    
    temp1 = (rows_educbin2 == 1);
    wvec_bins(i,2) = nanmean(sim_data(temp1,col_long.b)); % unconditional on CCT use
    vvec_bins(i,2) = nanmean(sim_data(temp1,col_long.x)); % unconditional on CCT use
    
    % educ bin 3
    temp1 = nansum(sim_data(rows_educbin3,col_long.CCT) == 1);
    temp0 = nansum(sim_data(rows_educbin3,col_long.CCT) == 0);
    yvec_bins(i,3) = temp1/(temp1 + temp0);
    
    temp1 = (rows_educbin3 == 1) & (sim_data(:,col_long.r) ~= 0);
    temp2 = sim_data(temp1,col_long.r);
    zvec_bins(i,3) = nanmean(temp2);
    
    temp1 = (rows_educbin3 == 1);
    wvec_bins(i,3) = nanmean(sim_data(temp1,col_long.b)); % unconditional on CCT use
    vvec_bins(i,3) = nanmean(sim_data(temp1,col_long.x)); % unconditional on CCT use
end

f = figure;
h = zeros(3,1); % initialize handles for 3 plots
h(1) = plot(xvec,yvec_bins(:,1),'k-o');
hold on;
h(2) = plot(xvec,yvec_bins(:,2),'k-d');
h(3) = plot(xvec,yvec_bins(:,3),'k-*');
hold off;
set(f,'visible','off');
xlabel('Child''s Age','Interpreter','LaTex')
ylabel('Fraction using CCT','Interpreter','LaTeX');
axis([0 17 0 1.1*max(max(yvec_bins+0.1))]);
h=legend('Father''s Educ. - High School or Less','Father''s Educ. - (Some) College','Father''s Educ. - Graduate','Location','southoutside');
set(h, 'interpreter','latex');
name = ['Fig_CCT_by_father_educ'];
saveas(f,name, 'png')
saveas(f,name, 'eps')

f = figure;
h = zeros(3,1); % initialize handles for 3 plots
h(1) = plot(xvec,zvec_bins(:,1),'k-o');
hold on;
h(2) = plot(xvec,zvec_bins(:,2),'k-d');
h(3) = plot(xvec,zvec_bins(:,3),'k-*');
hold off;
set(f,'visible','off');
xlabel('Child''s Age','Interpreter','LaTex')
ylabel('Average $r$ (if $CCT=1$)','Interpreter','LaTeX');
axis([0 17 0 1.1*max(max(zvec_bins+0.1))]);
h=legend('Father''s Educ. - High School or Less','Father''s Educ. - (Some) College','Father''s Educ. - Graduate','Location','southoutside');
set(h, 'interpreter','latex');
name = ['Fig_rpos_by_father_educ'];
saveas(f,name, 'png')
saveas(f,name, 'eps')

f = figure;
h = zeros(3,1); % initialize handles for 3 plots
h(1) = plot(xvec,wvec_bins(:,1),'k-o');
hold on;
h(2) = plot(xvec,wvec_bins(:,2),'k-d');
h(3) = plot(xvec,wvec_bins(:,3),'k-*');
hold off;
set(f,'visible','off');
xlabel('Child''s Age','Interpreter','LaTex')
ylabel('Average $b$','Interpreter','LaTeX');
axis([0 17 min(min(wvec_bins(:,:))) 1.1*max(max(wvec_bins(:,:)+0.1))]);
h=legend('Father''s Educ. - High School or Less','Father''s Educ. - (Some) College','Father''s Educ. - Graduate','Location','southoutside');
set(h, 'interpreter','latex');
name = ['Fig_b_by_father_educ'];
saveas(f,name, 'png')
saveas(f,name, 'eps')

f = figure;
h = zeros(3,1); % initialize handles for 3 plots
h(1) = plot(xvec,vvec_bins(:,1),'k-o');
hold on;
h(2) = plot(xvec,vvec_bins(:,2),'k-d');
h(3) = plot(xvec,vvec_bins(:,3),'k-*');
hold off;
set(f,'visible','off');
xlabel('Child''s Age','Interpreter','LaTex')
ylabel('Dollars per week','Interpreter','LaTeX');
axis([0 17 0 1.1*max(max(vvec_bins(:,:)+0.1))]);
h=legend('Father''s Educ. - High School or Less','Father''s Educ. - (Some) College','Father''s Educ. - Graduate','Location','southoutside');
set(h, 'interpreter','latex');
name = ['Fig_x_by_father_educ'];
saveas(f,name, 'png')
saveas(f,name, 'eps')

%% Graphs: Hourly wages and incomes as a function of PARENTAL age and education

% A) By parental age

xvec_l = [25:2:55]'; % lower bound ages
xvec_u = [26:2:56]'; % upper bound ages
yvec = NaN*ones(length(xvec_l),6);
for i = 1:size(xvec_l,1)
    age_lb = xvec_l(i);
    age_ub = xvec_u(i);
    rows_mom_data = data(:,col_long.mom_age)>=age_lb & data(:,col_long.mom_age)<=age_ub;
    rows_mom_sim = sim_data(:,col_long.mom_age)>=age_lb & sim_data(:,col_long.mom_age)<=age_ub;
    rows_dad_data = data(:,col_long.dad_age)>=age_lb & data(:,col_long.dad_age)<=age_ub;
    rows_dad_sim = sim_data(:,col_long.dad_age)>=age_lb & sim_data(:,col_long.dad_age)<=age_ub;
    
    % 1.a) mom hourly wage, data
    yvec(i,1) = nanmean(data(rows_mom_data,col_long.mom_wage));
    % 1.b) mom hourly wage, sim
    yvec(i,2) = nanmean(sim_data(rows_mom_sim,col_long.mom_wage));
    % 2.a) dad hourly wage, data
    yvec(i,3) = nanmean(data(rows_dad_data,col_long.dad_wage));
    % 2.b) dad hourly wage, sim
    yvec(i,4) = nanmean(sim_data(rows_dad_sim,col_long.dad_wage));
end

% 1. Mom hourly wage
f = figure;
plot(xvec_l, yvec(:,1),'k-',xvec_l, yvec(:,2),'k--');
set(f,'visible','off');
xlabel('Mother''s Age');
ylabel('Dollars per hour');
axis([min(xvec_l) max(xvec_l) 0 1.1*max(max(yvec(:,1:2)))]);
h = legend('Mother''s Wage $w_1$ - Data','Mother''s Wage $w_1$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_w1_age', 'png')
saveas(f,'Fig_w1_age', 'eps')

% 2. Dad hourly wage
f = figure;
plot(xvec_l, yvec(:,3),'k-',xvec_l, yvec(:,4),'k--');
set(f,'visible','off');
xlabel('Father''s Age');
ylabel('Dollars per hour');
axis([min(xvec_l) max(xvec_l) 0 1.1*max(max(yvec(:,3:4)))]);
h = legend('Father''s Wage $w_2$ - Data','Father''s Wage $w_2$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_w2_age', 'png')
saveas(f,'Fig_w2_age', 'eps')

% B) By parental education
xvec = [12:17]'; % possible levels of education
yvec = NaN*ones(length(xvec),6);
for i = 1:size(xvec,1)
    educ_i = xvec(i);
    rows_mom_data = data(:,col_long.mom_educ)==educ_i;
    rows_mom_sim = sim_data(:,col_long.mom_educ)==educ_i;
    rows_dad_data = data(:,col_long.dad_educ)==educ_i;
    rows_dad_sim = sim_data(:,col_long.dad_educ)==educ_i;
    
    % 1.a) mom hourly wage, data
    yvec(i,1) = nanmean(data(rows_mom_data,col_long.mom_wage));
    % 1.b) mom hourly wage, sim
    yvec(i,2) = nanmean(sim_data(rows_mom_sim,col_long.mom_wage));
    % 2.a) dad hourly wage, data
    yvec(i,3) = nanmean(data(rows_dad_data,col_long.dad_wage));
    % 2.b) dad hourly wage, sim
    yvec(i,4) = nanmean(sim_data(rows_dad_sim,col_long.dad_wage));
end

% 1. Mom hourly wage
f = figure;
plot(xvec, yvec(:,1),'k-',xvec, yvec(:,2),'k--');
set(f,'visible','off');
xlabel('Mother''s Years of Schooling');
ylabel('Dollars per hour');
axis([min(xvec) max(xvec) 0 1.1*max(max(yvec(:,1:2)))]);
h = legend('Mother''s Wage $w_1$ - Data','Mother''s Wage $w_1$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_w1_educ', 'png')
saveas(f,'Fig_w1_educ', 'eps')

% 2. Dad hourly wage
f = figure;
plot(xvec, yvec(:,3),'k-',xvec, yvec(:,4),'k--');
set(f,'visible','off');
xlabel('Father''s Years of Schooling');
ylabel('Dollars per hour');
axis([min(xvec) max(xvec) 0 1.1*max(max(yvec(:,3:4)))]);
h = legend('Father''s Wage $w_2$ - Data','Father''s Wage $w_2$ - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_w2_educ', 'png')
saveas(f,'Fig_w2_educ', 'eps')

%% Graphs: Non-Labor Income - plot average NLI by father's age and education, data vs. simulation

NLI_sim_data = fp.NLI_sim_data;

% A) By father's age
xvec_l = [25:2:55]'; % lower bound ages
xvec_u = [26:2:56]'; % upper bound ages
yvec = NaN*ones(length(xvec),2);
for i = 1:size(xvec_l,1)
    age_lb = xvec_l(i);
    age_ub = xvec_u(i);
    rows_dad_data = data(:,col_long.dad_age)>=age_lb & data(:,col_long.dad_age)<=age_ub;
    rows_dad_sim = NLI_sim_data(:,col_long.dad_age)>=age_lb & NLI_sim_data(:,col_long.dad_age)<=age_ub;
    % data
    yvec(i,1) = nanmean(data(rows_dad_data,col_long.nonlabor_income));
    % sim
    yvec(i,2) = nanmean(NLI_sim_data(rows_dad_sim,col_long.nonlabor_income));
end

f = figure;
plot(xvec_l, yvec(:,1),'k-',xvec_l, yvec(:,2),'k--');
set(f,'visible','off');
xlabel('Father''s Age');
ylabel('Dollars per week');
axis([min(xvec_l) max(xvec_l) 0 1.1*max(max(yvec(:,1:2)))]);
h = legend('Non-Labor Income - Data','Non-Labor Income - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_NLI_age', 'png')
saveas(f,'Fig_NLI_age', 'eps')

% B) By father's education
xvec = [12:17]'; % possible levels of education
yvec = NaN*ones(length(xvec),2);
for i = 1:size(xvec,1)
    educ_i = xvec(i);
    rows_dad_data = data(:,col_long.dad_educ)==educ_i;
    rows_dad_sim = NLI_sim_data(:,col_long.dad_educ)==educ_i;
    % data
    yvec(i,1) = nanmean(data(rows_dad_data,col_long.nonlabor_income));
    % sim
    yvec(i,2) = nanmean(NLI_sim_data(rows_dad_sim,col_long.nonlabor_income));
end

f = figure;
plot(xvec, yvec(:,1),'k-',xvec, yvec(:,2),'k--');
set(f,'visible','off');
xlabel('Father''s Years of Schooling');
ylabel('Dollars per week');
axis([min(xvec) max(xvec) 0 1.1*max(max(yvec(:,1:2)))]);
h = legend('Non-Labor Income - Data','Non-Labor Income - Simulated','Location','southoutside');
set(h, 'interpreter','latex');
saveas(f,'Fig_NLI_educ', 'png')
saveas(f,'Fig_NLI_educ', 'eps')
